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Abstract 

We propose a modification of the Hybrid Monte-Carlo method to sample 
equilibrium distributions of continuous field models. The method allows an effi- 
cient implementation of Fourier acceleration and is shown to reduce completely 
critical slowing down for the Gaussian model, i. e., z = 0. 
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The development of efficient numerical algorithms to study equilibrium properties 
of field-theoretical models near second order phase transitions is a very important 
subject in particle and statistical physics]!], [| f|, |^]. For that purpose, several 
methods such as Molecular Dynamics, Langevin and Monte-Carlo (MC) have been 
used. While the first two methods suffer from systematic step-size time discretiza- 
tion errors which may affect the computed mean values of observables, the only errors 
present in MC methods are of statistical origin and can be easily controlled, in prin- 
ciple, by varying the number of samplings. In practice, however, it turns out that for 
many problems of interest the number of configurations necessary to achieve a given 
small error is very large and grows as some power of the system size, thus requiring 
too much computer time. 

In its simplest form, MC introduces a stochastic dynamics which involves the 
proposal of a random new local field configuration plus an acceptance/rejection step. 
This method can be inefficient to implement in some systems. The Hybrid Monte 
Carlo (HMC) algorithm first proposed by Duane et al. |§, uses Molecular Dynamics 
to propose new configurations then requiring a single acceptance/rejection step every 
time the whole system is updated, which is a substantial improvement over other 
MC methods. The usual implementation of HMC relies on an appropriate numerical 
integration of the corresponding Hamiltonian dynamics of the system. 

In this paper, we present a generalization of HMC which is based upon the numer- 
ical integration of a non-hamiltonian dynamics which, however, conserves the energy 
of the system 0. Our method turns out to be closely related to the use of a time-step 
matrix in a Langevin integration scheme || and thus can be regarded as an exact 
version of this method in the sense of not being affected by step-size errors. We first 
define our generalized algorithm showing that it obeys correctly the detailed balance 
condition. As an application, we study the Gaussian model. The simplicity of the 
model enables us to carry a systematic analytic study of correlation times. We also 
show that the optimal implementation of our method for the Gaussian model is ac- 
tually a Fourier acceleration scheme which completely reduces critical slowing down 
(CSD). 

CSD theory|| tells us that near second-order phase transitions the correlation 
time, Tfi, of a measured observable O, increases as Tq ~ £ z , being £ the correlation 
length and z the dynamical critical exponent, ta can be defined as some measure of 
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the relaxation time of the correlation function of the observable O: 

< d(t)d(o) >-< d(t) >< d(o) > 
° < 6 2 (o) > - < 6(0) > 2 

In a numerical simulation of a d-dimensional system the computer time, Tq, required 
to measure the observable O at a given error behaves as T Q = t 6 L d ~ L d+C (2r 6 + 
1). The factor L d is always present as we simulate a system of size N = L d as 
a consequence of the increase of the number of degrees of freedom. The factor L c , 
present in HMC, is a consequence of the fact that simulations at a constant correlation 
length at bigger sizes may require an additional computer effort in order to keep 
the acceptance probability within reasonable limits. The last factor (2tq + 1) takes 
into account the number of effectively independent configurations produced by the 
algorithm. For finite systems close enough to the critical point, the correlation time 
Tq increases with system side as r Q ~ L z . For the local updating schemes such as 
heat-bath or Metropolis the exponent z being near 2 strongly demands on computer 
time (although the above defined exponent c is actually for these algorithms). For 
spin models the collective updating scheme of Swendsen and Wang [JlOj has proven 
quite successful in reducing the dynamical critical exponent and overcoming CSD. 
For continuous models a Multigrid Monte-Carlo Methodjll] was proposed which can 
also reduce CSD in certain cases. Also the time-step matrix Langevin method M 
can be helpful for some models if an appropriate matrix is chosen. The last two 
algorithms were shown to reduce the exponent z for the Gaussian model but it is 
not clear if some reducing can be achieved in interacting models such as the 4 
model. The standard implementation of HMC was not observed to ameliorate CSD 
for the 4 model |12| and for the Gaussian model it was shown that the algorithm 
reduces z = 1 if an appropriate tuning of the trajectory length is taken. For an 
application of the standard HMC to study the two dimensional XY model see |[14|| . 
Yet another algorithm also used for simulation of these models is Overrelaxation (see 
H and references therein) which was also shown to produce z = 1 for the Gaussian 



model|15|. The modified HMC we propose allows to reduce z = for the Gaussian 
model. 

We now describe the generalized HMC algorithm. Let us consider a system 
(0i, 02, • • • , 0aO = [0], of N = L d scalar variables, whose statistical equilibrium prop- 
erties are defined through the Gibbs factor exp(— H), by its Hamiltonian Ti[(j)]. The 
variables [0] are considered to be generalized coordinates and a set of conjugate mo- 
menta (p 1 ,p 2 , . . . ,Pn) = [p] associated with a kinetic energy 7i K = is 
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introduced. The variable pi can be in general a vector variable with D components, 
Pi = ip}, p^, ■ ■ ■ ,pf)- The total Hamiltonian is H = H + Hk- We propose the 
following dynamics: 

d6- D N 

^ = EE(-4 S M (2) 



dpj 
dt 



s=lj=l 
N 

Y^WiiFi ,s = l,...,D 



or, written in more compact vector notation: 

ji D 



ft = 5>v (3) 



8=1 

(A 11 ) 1 F ,s = l,...,D 



dt 

where the A s are some linear operators which can be represented as a matrix, and Fj 
represents the force as computed from the Hamiltonian — -^-H. The standard HMC 
substitutes the above dynamics by the Hamiltonian dynamics which can be obtained 
from the above set of equations by considering D = 1 and A equal to the identity 
operator. 

It is an essential property that can be easily verified that the proposed dynamics 
in equations fl3|) exactly conserves energy, i.e., dH/dt = 0. For the approximate 
integration of the previous equations of motion the "leap-frog" scheme can be used, 
introducing a discrete mapping [(p(t) , pit)] — ► [<p(t + St),p(t + St)] = G St ([<f>(t) , p(t)]) 
of phase space, dependent on the time step St chosen. The total energy, as a result 
of the time discretization used in the leap-frog scheme, is no longer conserved and its 
variation can be controlled by varying St. The leap-frog approximation reads: 

D /r.N 2 D 

0' = <P + StJ2AY + { -^-J2A s (A s ) T F([<f ) ]) (4) 

s=l A s=l 

V' s = P s + 6 {(A s ) T (Fm) + F(W})) 

We define yet another mapping obtained from n iterations of the previous mapping, 
i.e. Q = (G 5t ) n . The configuration obtained when one applies Q is then accepted or 
rejected in such a way that detailed balance is verified in order to sample the canonical 
distribution for the Hamiltonian TC. As in the standard HMC the momenta variables 
are refreshed after every acceptance/rejection step according to the Gaussian distri- 
bution of independent variables exp(— Hk)- The evolution given by Q (n leap-frog 
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steps) and the acceptance/rejection step constitute what is called 1 MC trial. De- 
tailed balance is obeyed if one requires G St to be time reversible and area preserving 
and if one accepts the new configuration with probability min[l, exp(— A7Y)], where 
AH = H{g{[<j),p})) - H([<f>,p]). The time reversibility G St {[(/)', -p']) = [<p, -p] can be 
easily verified to be obeyed by equations @. One can also proof that the area pre- 
serving property is verified by these equations. The above properties are satisfied for 
arbritrary matrices A s provided the associated mapping G 5t remains a one to one 
mapping in phase-space. We have defined a variety of HMC-type methods charac- 
terized by a particular choice of matrices A s . One can then choose the matrices A s 
that better suit a particular problem. We have shown before how a particular version 
of the above generalized HMC can be used to simulate conserved order parameter 



systems flEf]. 

We compare now our method with the one introduced in reference || based upon 
the numerical integration of a Langevin equation using a matrix time-step. This 
method is based upon the observation that the stationary probability distribution of 
the Langevin equation: 

9 -^l = + V2Ur) (5) 
or dpi 

is precisely exp(— H). Here &(t) are stochastic Gaussian random variables of mean 
zero and correlations < £i( r )£j( r/ ) >= ( % ( K r — r ')- The solution of the equation is 
approximated by: 



5H 



(6) 



j L "Vj 

Where is an arbitrary matrix and rjj is a Gaussian variable of mean zero and 
correlations < r]iT]j >= 5ij. This corresponds exactly to the one step leap-frog ap- 
proximation of the generalized HMC introduced above (equation (H) with D = 1) if 
we identify: (5t) 2 /2 = 5r and AA T = e. The main difference between the two meth- 
ods is the presence of an acceptance/rejection step in the generalized HMC absent in 
the numerical integration of the Langevin equation. In that sense, we can say that 
the generalized HMC method introduced in this paper makes exact (in the sense that 
averages are not biased by the choice of the time step) the numerical integration of 
the Langevin equation using a matrix time step introduced in reference ||. 

As an application we have considered the Gaussian model defined by the following 
Hamiltonian: 



N 



(7) 
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index i runs over the N = L 2 sites of a 2- dimensional square lattice, with periodic 
boundary conditions (a similar analysis can be carried out in any spatial dimension 
but we refer to the case d=2 for simplicity). Vl is the usual lattice discretized version 
of the gradient operator. This problem can be better analyzed in Fourier space. The 
total Hamiltonian H. in terms of the Fourier transform of fields and momenta space 
is: 



n 



N 

£ 

k=l 



^-\k\ 2 + \\h\ 



where ui k is given by uj\ = \i + 4(sin 2 (/c x /2) + sin 2 (/c 2/ /2)) and 4> k and p k stand for the 
fields and momenta variables in Fourier space. We choose the number of momenta 
variables associated to a given field equal to 1, D = 1. Suppose that we choose for 
the matrix A, generating the dynamics, a diagonal matrix in Fourier space. Then 
after n leap-frog steps, equation (ffl) implies that: 



u k (f) k (n8t) 


= M n k 


^fe0fc(O) 


p k {n5t) 




Pk(O) 



(9) 



for k 



N. Matrices are given by: 



Ml 



cos(n9 k ) sxn(n9 k )/ cos(9 k /2) 

cos(9 k /2) sm(n9 k ) cos(n9 k ) 



(10) 



where we have introduced 9 k = cos _1 (l — c|/2) and c k = A k u k 5t and A k denoting 
the diagonal elements of the matrix A in Fourier space. In this model the different 
modes evolve independently of each other and the evolution equations are linear (this 
is similar to the standard HMC with A = l |fL3[| ). For the stability of the leap-frog 
integration, the eigenvalues of M k should lie on the unit circle of the complex plane 
which happens to be the case if c k is between and 2. 

By using the evolution equations together with the assumption that the field vari- 
ables 4> k (0) are in thermal equilibrium and, therefore, follow the distribution exp(— Tl), 
one can compute the equilibrium average discretization error as: 



N 



< H(n5t) - H(0) >=< AH >= ^ 



k=\ 



32-8 



sin 2 (n^ 



(11) 



In reference |T^] it was shown that, to a good approximation, the average ac- 
ceptance probability < pa > is related to the average discretization error < AH > 
by: 

< PA >= erfc(i v // < Aft >) (12) 
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We now turn to the question of the optimal choice for the matrix A. From 
equation ([D it is immediately seen that if we choose the matrix A such that A k = 
1/uJk the iteration equations get independent of the mass fi and all the modes are 
equally updated. This is, in effect, an exact implementation of the method of Fourier 
acceleration. This choice of the matrix clearly reduces completely CSD (z = 0) in the 
sense that correlation times are independent of the mass even when the mass goes 
to zero and the model becomes critical. The standard HMC corresponds here to the 
choice Ak = 1, independent of k [|13| . 

Let us compute the computational effort needed to achieve a given statistical error. 
In order to make further analytical calculations, it is convenient to introduce a set of 
random variables cr m which take the value 1 or if the configuration proposed after 
m MC trials has been accepted or not, respectively. Using this variable we can write 
an expression for the field variable after m MC trials, 4> k {m n8t), as 



<^fc0fc(^ nSt) = o Tl 



cos{n8 k )uj k </) k ((m - 1) nSt) + ^77r%r Pfc((ra - 1) n5t) 

cos(0 k /2) 



+ 



(1 - a m )u k (j) k ({m - 1) nSt) (13) 

The momenta p k (m nSt) are the independent random variables following a Gaussian 
distribution which are drawn after the m-th acceptance/rejection step. This equation 
can be iterated to obtain (f) k (m nSt) in terms of <fi k (0) and all the momenta generated 
during the evolution. 

The variables a m are Bernoulli variables with probability of being equal to 1 equal 
to the acceptance probability, min[l, exp(— AH(m)}. This probability depends on the 
total change in energy at the m th MC trial, ATi.(m) which is a function of the initial 
field configuration and of the momenta generated at the j th MC trial and variables <jj 
such that j < m. To proceed we make the approximation that the a variables are all 
independently distributed variables with the probability of being equal to one equal 
to the average acceptance probability p =< Pa >, given by flT2|) . This means that we 
consider the probability to accept or reject the whole configuration at a given step 
to be independent of the previous "time-history" of the system. The approximation 
is reasonably good as we will see later. Within this approximation, the correlation 
function for the magnetization is: 

2 77,00, 



Cm{iti) = Cjvf(m n8t) 



2p sin z . 
F v 2 



(14) 



The obtained correlation function is exponential with a correlation time equal to 
t m = — l/log(|C , Af(l)|) in units of MC trials. The relaxation time of other modes is 
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obtained replacing in the previous equation 6 by 6k- For the optimal matrix, 6k = 6 
is independent of k and so all the modes relax in the same way. 

In figure 1 we compare Cm(1) given by the above correlation function for the 
magnetization with simulation results for the case L = 32 and n = 4 as a function 
of St. The excellent agreement between the analytical expression and the simulation 
results shows that the previous approximation works extremely well for the calculation 
of the correlation function for the magnetization. The agreement between simulation 
and this approximation actually improves with increasing system size. We note the 
following features (see figure 1): the value of tm as a function of St for a given n 
shows a minimum for small values of St and then n — 1 zeros for n odd and n zeros for 
n even. This minimum disappears for n larger than a given value. We discuss now 
how the parameters n and St should be optimized in order to minimize the computer 
effort for the magnetization tu = {2t~m + l)n. 

In principle, since there is always a value of St that yields tm = for n = 2, these 
are obviously the optimal choices. The corresponding computational effort tu = 2 is 
independent of the system size corresponding to an exponent c = 0. However, a closer 
look shows that as the system size grows the precise St value at which the correlation 
time is zero is very difficult to locate in the sense that a slight error in the chosen value 
puts the system in a region of high (negative) correlations (see figure 1). Furthermore 
the corresponding correlation times of the energy would not be optimized, increasing 
enormously as the system size grows. This becomes apparent when one discusses the 
correlation function of the energy. The obtained approximation for the correlation 
function of the energy is: 

C H {m) = C H {m nSt) = [1 - psm 2 {n6)] m (15) 

The correlation time is = — 1/ log(C'^(l)). In figure 2 we have compared CV^(l) 
given by the above correlation function to simulation results. The agreement here, 
although still reasonable, is not as good as it was in the case of the magnetization. 
However, one can still obtain precise quantitative conclusions from this figure. The 
value of Crt(l) has no zeros as a function of St for a given n, in contrast to what 
happens with C'm(I)- It still shows, however, an absolute minimum at a given St 
which approaches zero as n is increased. We have used this minimum in the correlation 
function of the energy as a function of St to find the optimal n for a given system 
size that minimizes the computer effort (2r n + l)n. Although we have been unable 
to obtain an analytical expression for the optimal values for n and St, a numerical 
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study of equations (|ITD, (|l"2j) and ( |i5|) allows us to conclude that the optimal value 
for n increases with system size as L 1//2 whereas the optimal St behaves as L~ x l 2 . 
The obtained optimal values of n were such that Cm (1) still has the local minimum 
as a function of St mentioned above and which is near the absolute minimum of the 
energy. For large L, the corresponding correlation times Tm and turn out to be 
L-independent and are given by tm = 2.5, Tu = 1.5 approximately. The optimal 
acceptance probability p ~ 0.67 and the product nSt ~ 1 are also independent of the 
system size. 

In summary, the computational effort both for the magnetization and the energy 
behaves as t M , t n ~ L 1//2 , corresponding to an exponent c = 1/2. This can be 
understood in the following way: In order to keep the acceptance probability constant 
as we increase the system size, St has to be varied as L~ d ^. Thus one needs to increase 
n as L d / 4 in order that the product nSt appearing on the correlation function remains 
also unchanged as the system size grows. The same picture was also seen to apply to 



a study of the <f) A model performed with the standard HMC method ^] . One should 
also note that an explicit implementation of the method of Fourier acceleration needs 
the calculation of Fourier transforms that involve an additional computer effort of 
order logL. 

In conclusion, we have proposed a generalized version of HMC that is an ex- 
act implementation of time-step matrix Langevin methods. We have considered the 
Gaussian model as a test-case of the algorithm. The optimal matrix A is a diagonal 
matrix in Fourier space with diagonal elements = 1/uk (Fourier acceleration). 
This matrix reduces completely CSD in the sense that the correlation times are mass 
independent, z — 0. We have proposed an approximation for the calculation of the 
correlation functions of energy and magnetization that improves as the system size 
gets larger. Using this approximation we have discussed the optimization of the pa- 
rameters of the algorithm. The optimal computer effort grows as L d/4 due to the 
decrease of the acceptance probability with the system size and the need to keep the 
trajectory length, nSt, constant at different system sizes. 
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Figures Captions 



Figure. 1.- Comparison of Cm(1) as obtained from our analytical approximation 
(continuous line) from equation (0), for a system of size L = 32 and n = 4 with 
simulation results (rhombi). 

Figure. 2.- Comparison of Ch(1) (continuous line) from equation (|1|), with 
simulation results (rhombi) for the same system size and number of leap-frog steps 
used in figure 1. 
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